Stability of the puncture method with a generalized BSSN formulation 
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The puncture method for dealing with black holes in the numerical simulation of vacuum space- 
times is remarkably successful when combined with the BSSN formulation of the Einstein equations. 
We examine a generalized class of formulations modelled along the lines of the Laguna-Shoemaker 
system and including BSSN as a special case. The formulation is a two parameter generalization of 
the choice of variables used in standard BSSN evolutions. Numerical stability of the standard finite 
difference methods is proven for the formulation in the linear regime around flat space, a special case 
of which is the numerical stability of BSSN. Numerical evolutions are presented and compared with 
a standard BSSN implementation. Surprisingly, a significant portion of the parameter space yields 
(long-term) stable simulations, including the standard BSSN formulation as a special case. Fur- 
thermore, non-standard parameter choices typically result in smoother behaviour of the evolution 
variables close to the puncture. 

PACS numbers: 04.20.Cv, 04.25.D-, 04.25.dg 



I. INTRODUCTION 

Accelerated bodies generate gravitational waves 
(GWs) in analogy to the emission of electromagnetic 
waves by accelerated charges. The first direct detection 
of GWs, expected in the course of the next decade, will 
not only provide us with the first strong field tests of 
Einstein's general relativity but also open up an entire 
new window to the universe. The strongest source of 
GWs are compact binary systems involving neutron stars 
and black holes (BHs). Such compact objects have been 
known for a long time to represent the natural end prod- 
uct of stellar evolution. For instance, stellar-mass BHs 
are suspected to be the compact members in X-ray bi- 
naries In addition there is now strong observational 
evidence for the existence of supermassive BHs (SMBHs) 
at the center of many if not all galaxies 0, H| . Astrophys- 
ical observations in recent decades have thus promoted 
BHs from the status of a mathematical curiosity to that 
of a key player in many astrophysical processes. 

While GW emission from compact objects has been 
theoretically predicted for quite a while, the waves' weak 
interaction makes their direct observation a daunting 
task, possible only by using modern high precision tech- 
nology. In particular, there exists now an international 
network of ground-based laser interferometers (LIGO 
0, H, GEO600 H, VIRGO @ and TAMA % op- 
erating at or near design sensitivity. A space-borne in- 
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terferometer called LISA [T(| is scheduled for launch in 
about one decade to supplement such observations with 
exceptional accuracy in a lower frequency band. Still, the 
understanding of the radiated wave patterns is crucial for 
the first detection of GWs as well as for the interpreta- 
tion of the measured signal. Eventually, the community 
will be able to gain information about characteristic pa- 
rameters of the BH system observed via GWs such as the 
mass ratio and spins. 

The modeling of these binary sources of GWs currently 
employs a variety of techniques. The inspiraling phase of 
a binary black-hole (BBH) prior to merger as well as the 
ringdown phase after the merger can be modeled accu- 
rately by the approximate post-Newtonian and per- 
turbation methods [l2[ , respectively. In order to simulate 
the late inspiral and merger of a BBH, however, numer- 
ical methods are required to solve the fully non-linear 
Einstein equations. A numerical treatment requires us to 
cast the Einstein equations into the form of a time evo- 
lution system. This is most commonly done by using the 
canonical Arnowitt-Deser-Misner "3+1" decomposition 
[l3j as further developed by York [l4| ; the 4-dimensional 
spacetime is decomposed into a family of 3-dimensional 
hypersurfaces labeled by a time coordinate. The geom- 
etry of spacetime is determined by the induced 3-metric 
jij on the hypersurfaces and their extrinsic curvature 
Kij, which describes their embedding. The coordinates 
are described by the lapse function a and the shift vec- 
tor (3 l . These gauge functions represent the coordinate 
freedom of general relativity (GR). 

For a long time, numerical methods based on this ap- 
proach faced a variety of problems including the specific 
formulation of the evolution equations, suitable coordi- 
nate choices and the treatment of singularities inherent in 
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the spacetimes. The year 2005 brought about the even- 
tual breakthrough, the first complete simulations of a 
BBH coalescence |15l - lf7| . The ensuing years have pro- 
duced a wealth of results on BBH inspiral pertaining to 
BH recoil, spin precession and GW data analysis to name 
but a few (see [18l - l22l ] for recent reviews). 

The current generation of successful numerical codes 
can be divided into two categories. The first class uses 
the so-called Generalized Harmonic Gauge (GHG) formu- 
lation employed in Preterms' original breakthrough. The 
second type of codes is commonly referred to as Moving 
Puncture codes, the method underlying the simulations 
of the Goddard and Brownsville groups. In spite of the 
remarkable robustness of both methods, it is fair to say 
that our understanding of why these techniques work so 
well is limited. The Moving Puncture method in par- 
ticular has proven robust in even the most demanding 
simulations of BBHs involving nearly critical spins and 
velocities close to the speed of light [23l426l |. Previous 
investigations of this method have concentrated on the 
structure near the singularity and the impact of gauge 
conditions l27l-[3ll. 



The purpose of the present work is to shed addi- 
tional light on which ingredients of the Moving Punc- 
ture method make this technique so successful. The par- 
ticular focus of our study is on the underlying formula- 
tion of the Einstein equations, the Baumgarte-Shapiro- 
Shibata-Nakamura (BSSN) formulation [32|, |33| as well 
as a modified version thereof modeled along the lines of 
the alternative evolution system proposed by Laguna and 
Shoemaker (LaSh) in 2002 [34| . Such a study is beyond 
purely academic interest. While the currently employed 
techniques appear to work well for 3+1 dimensional sim- 
ulations in the framework of Einstein's general relativity, 
there is strong motivation to push numerical relativity 
further. A main target of gravitational wave observations 
is the testing of GR versus alternative theories of grav- 
ity (see 35 1 for an overview, [36j | for solutions of rotating 



holes in Chern-Simons modified gravity and [37J for hy- 
perbolicity studies of scalar tensor theories of gravity). A 
further application of numerical relativity, in the context 
of high energy p hysics, as motivated by so-called TeV 
gravity scenarios |38l442| . or by the (conjectured) Anti- 
de Sitter/ Conformal Field theory (AdS/CFT) correspon- 



dence [43W45I . will be the simulation of BHs in higher 
dimensions or non-asymptotically flat spacetimes 

(see, e.g., [561 ] for a recent approach). An improved un- 
derstanding of the success of the 3+1 GR techniques will 
be crucial in extending numerical relativity successfully 
along these lines of future research. From a more prac- 
tical point of view, alternative schemes might simply be 
more efficient and result in reduced computational re- 
quirements. Unfortunately, we will see further below that 
the LaSh system does not result in faster simulations. 

This paper is structured as follows. The formulation 
of the LaSh evolution scheme is presented in Sec. [TTJ 
In Sec. IIIII well-posedness and numerical stability of the 
LaSh system are studied. The LaSh formulation, imple- 



mented as an extension to the Lean code [57[, is tested 
numerically with head-on collision and inspiraling BH bi- 
naries. The numerical results are presented in Sec. IIVI 
Finally Sec. fVl contains our conclusions. 



II. THE LASH FORMULATION 
A. The ADM equations 

Both the BSSN and LaSh systems are typically pre- 
sented as a simultaneous conformal decom pos ition and 
readjustment of the ADM equations (33l, [58l,l59j . For our 
purposes such a presentation will not suffice. Instead, 
the addition of definition-differential constraints which 
alters the characteristic structure of the system and guar- 
antees well-posedness and the conformal decomposition 
that changes to a convenient form of the evolved variables 
are considered separately. 

In any case one must first introduce the ADM system, 
which has the evolution part 

dtjij = -2aKij + Cpjij, (1) 
d t Kij = -D i D j a + a{R ij -2K ik Kf+K ij K] 

+C P K 13 . (2) 

and the physical Hamiltonian and momentum constraints 



H = R + K 2 = 0, 

Mi = DjKj - D t K = 0, 



where 



p& yk 

-'H? ~~ 1 ij,k 1 kj,i 



rk Tit 
kl L ij 



1 il l kj- 



(3) 
(4) 



(5) 



When closed with some gauge choice, the ADM system is 
typically only weakly hyperbolic and thus does not admit 
a well-posed Cauchy problem. The BSSN formulation is 
one of many modifications to the ADM system that can 
yield a strongly (or even symmetric) hyperbolic problem 
when coupled to some gauge [6C| . 



B. BSSN Constraint addition 

Definition- differential constraint: We define the dif- 
ferential constraint 



G* 



(G) 



Below it will be seen that this choice naturally makes fi 
coincide with the relevant BSSN variable. 



3 



Constraint addition: The ADM equations are ad- 
justed to 



dtla = ADM, 



d t K i3 = ADM + aG {lJ) 



(7) 
(8) 



+ 2aM t - 2aG j K\l + C p G l 
+ ~ fij G k d k p i -^GiBjp* , 



ADM 



(9) 



where tf denotes the trace-free part. The principal part, 
i.e., highest derivatives of variables added correspond ex- 
actly to those added in the Nagy-Ortiz-Reula (NOR) for- 
mulation [60 | (w ith a = b = 1, c = d = — 1/3 in the 
notation of |61(). It is for this reason that Gundlach 
and Martm-Garci'a were able to identify the two systems 
when analyzing the principal part (6l|V 



evolve according to 

dtnu = -2a X 3nK/2 A l3 + p% jtk + 2%^ 

- (16) 

dtX = /3 l X, 4 + \x{aX inK K - (17) 

d t A\j = x ~ 3nK/2 [-If i D j a + atfjf 

+ (1 - n K ) X 3nK/2 aKA^ + p k A i j , k 

- A k 3 P] k + A* k (3 k + ukA^/3% , (18) 
d t K = - X - 3/2nK D t D*a + (5 k K. k + n K K/3 k k 

+ X 3nK/2 a(A^A l3 + (1 - 3n K )K 2 /3), (19) 
d t r = -2 X 3nK/2 A tj a^ + 2a{ X 3riK/2 T) k A lk 

- \ X 3nK/2 A^ HX)J ~ \l l] {x 3nK/2 K), 3 ) 



C. Conformal decomposition and densitization 



Conformal variables and algebraic constraints: The 
LaSh system [34] takes as its evolved variables 



7ij 


= 7 3 lij, 


(10) 


X 


1 

= 7 3 , 


(11) 


K 


= X~ inK K, 


(12) 


A 1 ■ 

A 3 


= X -> n *(K i j -5 i j K/3), 


(13) 


F 




(14) 



The key difference between LaSh and BSSN is that inside 
LaSh the trace and tracefree parts of the extrinsic curva- 
ture are densitized. Notice, that we recover the standard 
BSSN equations for vanishing densitization parameter. 
Note that the definition of Gi gives 



Gi = fi- 7 jfc (lij,k ~ g7jfe,i 



= 7^ - J jk %, k . 



(15) 



Evolution equations and constraints: Taking a time 
derivative of the definitions, substituting the evolution 
equations and rewriting in terms of the evolved variables 
gives the LaSh equations - up to the algebraic constraints 



D 



ln(det 7 ) = 0, S = %A l j} = 0, T = f J A 



0, 



which are assumed to be satisfied exactly. The unknowns 



where ( J, denotes the definition of those terms rather 
than the evolved variable and Rij is partially rewritten 
in terms of T 1 , 



R i3 ~ R ij + Rij > 



1 ~ / ~ \ fc ~ 

~,ltn~, i ~ j^k 



R *3 = ~ ~1 m iij,lm + 7fc(i| r | ,j) + ( F ) , r (u)fc 



I 1 ' 



( 2 rf(if :) -)fc m + f k m tkij 



(21) 



(22) 



(23) 



Di denotes the covariant derivative compatible with the 
conformal metric. The physical constraints are rewritten 



H = R— x 3nK {A j k A k j - \k 2 ) 

Mi = Ah,i - \k,i - —k x ,i - timf ™ 

-i-(l-n K )I m iX ,m = 0. 
2x 



0, 



(24) 



(25) 



The differential constraints are given by Eq. [15] and alge- 
braic constraints by 



SsTipi'j-j =0, 
D = ln(det 7 ) = 0. 



T = fiAij = , 



(26) 
(27) 



Numerical relativity codes use a technique called con- 
straint projection to enforce the algebraic constraints. 
When operations are performed which may violate the 
D, S and T constraints they are enforced explicitly. It is 
for this reason that we need not worry about the algebraic 
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constraints in the construction of (|16fl20|l ; the continuum 
system they represent is identical to that of (JTJIHl). BSSN 
evolves Ay , so does not have the symmetry constraint S. 



D. Gauge conditions 

The successful evolution of binary BH systems has 
been possible with the now standard 1+log variant of 
the Bona-Masso slicing condition, 



d t a = -2aK + /3 l d l a. 



(28) 



Stationary data for this gauge has been studied in [27 
HH, l62| . In our numerical evolutions the T-driver shift 
condition 

9 t p i = HsB i + ZiP>d j i , (29) 
8 t B l = d t f l - bpdjP - nB l + ti^d 3 B\ (30) 

is used with {us, £i, £2, v) = (1, 0, Q, 1) unless otherwise 
stated. We refer to the combination of the 1+log lapse 
and T-driver shift as "puncture gauge". Conditions in 
which the lapse and shift are promoted to the status of 
evolved variables are often called live gauge conditions. 
In analytic studies however it is common to consider a 
fixed, densitized lapse 



Q = i 



"Q 



X 



(31) 



and fixed shift. In contrast to BSSN the LaSh system 
takes the densitized lapse Q as a dynamical variable. It 
is evolved according to the 1+log condition (|28|) rewritten 
in terms of Q. The original LaSh system [3J] is modified 
here by the consideration of different densitization pa- 
rameters riK and uq for the extrinsic curvature and the 
lapse. When comparing the computational cost of BSSN 
and LaSh in section lTVl we additionally evolve LaSh in the 
downstairs form of the conformal extrinsic curvature Aij . 
Whereas in 63] the focus was on changing the partial dif- 
ferential equation properties of the formulation and hold- 
ing the variables fixed, here we consider the effect of a 
change of variables alone. 



the theoretical background. Next, in Sec. MI Bl we use 
characteristic variables to demonstrate well-posedness for 
the continuum system. The analysis is then extended to 
the semi-discrete case, and follows closely the method of 
[fril lo5| . Finally, we deal with the algebraic constraints 
in the fully-discrete system by demonstrating that the 
natural semi-discrete limit of the standard implementa- 
tion (with constraint projection) is given by the systems 
considered in Sec MI Bl 



A. Theoretical background 

Continuum system: The linear, constant coefficient, 
first order in time and second order in space time evolu- 
tion problem 



d t u = P[d x ]u, u(t = 0, x) = f(x) 



(32) 



is called well-posed with respect to a norm || • || if for 
every smooth, periodic f{x) there exists a unique smooth 
spatially periodic solution and there are constants C, K 
such that for t > 



Mi,-)||<^e Ct || W (0 



(33) 



A hermitian matrix H(ui) is called a symmetrizer of the 
system if the energy u* Hu is conserved by the principal 
part of the Fourier transformed system, with 



w a 
/ / ' 



(34) 



for some K > constant, for every frequency uj in Fourier 
space (u denotes the Fourier transformed function.) We 
say that the Hermitian matrices A, B satisfy the inequal- 
ity A < B if Ay < By for every y. Well-posedness is 
equivalent to the existence of a symmetrizer (6(| , which 
is in turn equivalent to the existence of a complete set of 
characteristic variables. 

Discrete system: We introduce a grid 



x j = fepfe)^) = (jihjzhjah), 



(35) 



III. WELL-POSEDNESS AND NUMERICAL 
STABILITY 

The well-posedness of the LaSh system with either the 
puncture gauge or a fixed densitized lapse and shift (as- 
suming that the D,S and T constraints are satisfied) was 
previously studied in [58|, [6l[, so verifying these prop- 
erties for the system linearized around flat space is in 
its own right uninteresting. However we wish to demon- 
strate the numerical stability of the LaSh system around 
flat-space. The approach for the semi-discrete scheme 
is analogous to that for the continuum system, so we 
first tackle that problem. In Sec. MI Al we briefly recap 



with ji = 0, . . . N r —1 and h — 4£ is the spatial resolution. 
We denote 

D+Vi = ^-{vi+i - Vi), D_Vi=\(vi- t>i_i), (36) 
h h 



1 . 

1 = 2/i i+1 ~ Vj - 1 >- 



(37) 



The standard second order accurate discretization is writ- 
ten 

dt -> D Qt , (38) 

d idj ^ D% = { ° oiD J>i (39) 
J w \ D +i D-i 1 = 2 y ' 
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For brevity, we do not consider higher order discretiza- 
tion. Fourier transforming reveals 



m = 



i = i 



(40) 
(41) 



We abbreviate Sj = sin£j and = sin£j/2, and write 

2tt 



£ r = US r h = — 7T - 



N 



-in 
TV' 



(42) 



where uj r = — ^ + 1, and r = 1,2,3. The time 
step k is related to the spatial resolution h through the 
Courant factor k — X c h. Finally, we use the notation 



or 



-r] ij D m D 



0j) 



n 2 



-V ij D. 



(2) 
ij ; 



(43) 



where if 3 is just the identity matrix, so that J7 2 = 
53j=i |-D+i| 2 - The results for numerical stability with 
a polynomial method of lines time-integrator are anal- 
ogous to the result at the continuum: if there exists 
a hermitian H(£) for every grid frequency £ such that 
the energy u* Hit, is conserved by the Fourier-transformed 
semi-discrete principal system and satisfies 



a system guarantees the existence of a complete set of 
eigenvectors, which in turn guarantees well-posedness in 
some norm. For the discrete system, a similar approach 
is not possible with the standard discretization because 
the various blocks of the system remain coupled. This 
complication is caused by the fact that under the stan- 
dard discretization the second derivative is not equivalent 
to a repeated application of the first derivative. We will 
see in the following sections that this forces us to con- 
sider significantly larger matrices, and that for the main 
case of interest, the stability of the LaSh formulation 
with puncture gauge, the calculation is impractical. One 
may also consider the numerical stability of systems with 
the Dq discretization, in which second derivatives are ap- 
proximated by repeated application of the centered dif- 
ference operator Dq. In this case it is possible to make a 
2 + 1 decomposition of the semi-discrete system. Unfor- 
tunately, the discretization suffers from the problem that 
the highest frequency mode on the grid is not captured 
by the scheme. For the Fourier transformed system this 
property implies that the transformed spatial derivatives 
vanish, which typically prevents one from building an 
estimate on the highest frequency mode. Although arti- 
ficial dissipation may restore stability, we do not consider 
the Dq discretization further. 



K~ L In <H<KIr, 



In = 



n 2 
o 



(44) 



with K as above, then it is possible to construct a discrete 
symmetrizer for the semi-discrete problem without lower 
order terms. If the spectral radius of the product of the 
time-step and the semi-discrete symbol is bounded by 
a value that depends on the time-integrator, then the 
system is stable with respect to the norm 



MIm>+ = IWIh + ll«/llh + Z)l 



l^+^lll 



(45) 



where the subscript distinguishes between variables that 
appear as second derivatives in the continuum system. 
The estimate 



iAt| 



\h,D- 



< Ke 



CnAt 



iKlk 



D+, 



(46) 



then holds. For details we refer the reader to [64j |. 

Discussion: A straightforward way to construct char- 
acteristic variables for the continuum system is to per- 
form a 2 + 1 decomposition. One then ends up with de- 
coupled scalar, vector and tensor sectors which are hope- 
fully straightforward to diagonalize. Diagonalisability of 



B. Continuum system 

Fixed densitized lapse and shift: We begin by lineariz- 
ing the LaSh system around flat-space. Following [64ll66j 
we Fourier transform in space, and make a pseudo- 
differential reduction to first order. Spatial derivatives 
transform according to di — > tWi. The system has a 
complete set of characteristic variables with character- 
istic speeds (0, ±w, i^ngw). A conserved quantity for 
the system may be trivially constructed from the char- 
acteristic variables. It is straightforward but tedious to 
demonstrate that the conserved quantity is equivalent to 
the norm 



Mlfd 



lb* 



\Kt. 



\\h 



Ell7«, 



(47) 



k=l 



Puncture gauge: For simplicity we consider the time- 
integrated T-drivcr shift condition 



m - ft 



(48) 



The transformed vector of evolved variables is u 
(%,&, fk, Ki m , fin)- The principal symbol is 



6 



5*' 



tf 



eft ltf 



-2S\5f 
-2rf m 






2ztuw(j5^ 


: (5," + &"u> k ) 





(49) 



We denote Wi — wwi, with w — \w\. Here "trace- free" 
denotes that the trace is removed in downstairs indices. 



The characteristic variables can be constructed from the 
matrix 



J 










—iw 




y ±iw[8 P 5 J q + ia) p d) 9 77 ij ] tf -^iw[u) p &, 



itf 



— &(i<5j) + LU i UJ j UJ k 







2 lm 
3 





I cm 
pOq 



| tf ±2iw[w„w a l tf cj n 



/ 

(50) 



through U c = T~ lr u. The characterist ic sp eeds corre- 
sponding to each row are (0, ±y2w, ±^2/3w, ±u>, ±w). 
The conserved quantity is given by 



Parseval's relation implies equivalence with 



Il7d| 2 



\K, 



ei 



1 

£3 



| + - + w 2 |d| 2 



aei 



/ e 3 + C6 + 



2e 7 



4a 



+ / 



28 
~9~ 



256 



f I n + 2 +d+ ^7 f ' W H 



1 4 

£4 e 6 

..2 I ' |2 



8 

3es 



ha 



(51) 

2 



(2rf + 4e+16/)w 2 |A| 2 + (3c + 2d+^/) \K l3 \ 2 



436 

+ ^a(3 + e 2 ) + 3rf + 4e + / ^14 + £4 + e 6 + ^) ) |/i| 2 - 

By choosing a = kj, c = d = / = 1, e = 26 and £1 = 2, 
£2 = £4 = £6 = £8 = I, e 3 = 4, £5 = 16, £7 = 8 we obtain 



■ft^lMlL < e c < K\\ii\ 12 



Ipg - tl ii"'iipg (52) 

where = 125, and have demonstrated that the con- 
served quantity is equivalent to the norm 



ipg 



|7«H a +«' 2 ||&l| 2 +«' 2 ||Al| 2 + ||^«|| a + ||/i|| a . 

(53) 



2 1 M„ 1 1 2 

+ \ \ a ,k\\ 



H,k\ 



fe=l 



in physical space. 



(54) 



C. Discrete system 



Fixed densitized lapse and shift: We now consider 
the semi-discrete system with fixed densitized lapse and 
shift. As in the continuum case we linearize around flat- 
space. The difference operators transform as described 
in Sec. IIII Al We consider only the case tiq = 1. Wc 
define 



= .u - -D 0i r. 



(55) 



Decomposing the system into trace, off-diagonal, and di- 
agonal terms adjusted by the weighting tfju = ju and 
tfTi = Ti various sectors of the system decouple. In the 
following jij explicitly means i ^ j. The principal sym- 



7 



bol is 





^ 2 







f2 2 1 







if) 2 1 



The characteristic variables are 



K ± — T, 

2 ' 




if) z 



tf 



(56) 
(57) 

(58) 
(59) 

(60) 
(61) 

(62) 
(63) 



and have speeds (0, ±ifl, ±zf2). The system has a 
pseudo-discrete reduction to first order that admits a 
symmetrizer for every grid-frequency. One must treat the 
lowest frequency separately, but in that case the prin- 
cipal symbol vanishes and so admits the identity as a 
symmetrizer. By the equivalence of norms in finite di- 
mensional vector spaces we then have numerical stability 
in the pseudo-discrete norm 



Mhjd = ^\\%\\h + + WMi + \\%\\t (64) 

provided that the von-Neumann condition given by 



2X2 



(65) 



where Co = 2 and Co = v8 for iterated Crank-Nicholson 
or fourth-order Runge-Kutta and 2%2 = Qh, is satisfied. 



Parseval's relation guarantees equivalence with the dis- 
crete norm 



W\\t D+ j d =\\l l3 \\l + \\K l3 \\l + \\M\l+ 



E 

k=l 



\D + klr 3 \\l 



(66) 



in physical space. 

Puncture gauge: The principal symbol is a 19x19 ma- 
trix which contains several parameters. We are able to 
compute characteristic speeds for the system, but they 
are complicated, so we do not display them here. The 
lowest freqency mode is again trivial to analyze. In that 
case the pseudo-discrete reduction to first order has a 
vanishing principal symbol, and thus admits the iden- 
tity as a symmetrizer. For the non-maximal modes the 
principal symbol is complicated. As previously stated, 
performing a 2 + 1 decomposition on the semi-discrete 
symbol is not helpful, since the various sectors of the sys- 
tem remain coupled. We were therefore unable to find 
the eigenvectors of the matrix. We considered various 
subsectors of the full system. Since the (a,K) subsector 
is exactly the second order in space wave equation, it is 
trivial to demonstrate numerical stability. We considered 
also the subsector (fi,j3j) with the other variables frozen, 
and find characteristic speeds and variables. Once the 
two blocks are coupled to give the (aJi,K,(3j) subsec- 
tor, we did not manage to compute eigenvectors in finite 
time. It is possible to find a complete set of characteristic 
variables for the highest frequency grid-mode, since the 
principal symbol in that case takes a simpler form. In 
Sec. lIVI we present robust stability tests of the numerical 
implementation, which provide evidence that the system 
is formally numerically stable. 

Puncture gauge with non-standard spatial discretiza- 
tion: We are able to find a complete set of characteris- 
tic variables for a slightly altered spatial discretization. 
If one insists on using the Dq operator for the divergence 
terms didj/3 J in the evolution of fi the principal symbol 
becomes, with u =(%,&, fk,K im , /?„), 



P£ = 



\ 












-26lS™ 


2Ao^") 











-2r/ tm 
















-fl 2 <5£ + lD%D ok 


3 ^Im J 


-A (2) 

























(67) 



As before, one has to consider the lowest and high- ferent form. For the lowest frequency mode the prin- 
est frequency grid modes seperately, because in those cipal symbol of the pseudo-discrete reduction to first 
cases the principal symbol of the system takes a dif- order again vanishes, and can be dealt with as before. 
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For the sub-maximal frequencies the characteristic vari- 
ables can be constructed and have characteristic speeds 
(0, ±t^/2Q, ±z^/(3fi 2 + w 2 )/3, ±ifl, ±itt). The 0-speed 
characteristic variable is 

U = Doif + X -(u? - n 2 )a + i(3r! 2 + u?)t. (68) 
The lapse characteristic variable is 

U ±V z = V2K ± iVLa, (69) 
The longitudinal shift characteristic variable is 



4/ 



4%/3 



K± 



V3 



rA)if 



V3ft 2 +w 2 (w 2 - 3ft 2 

and the transverse shift modes are 

U i± i =n(D 0i D% +Lo 2 8 k l )p k 
±(D 0l D k +u J 2 6 k t )f k . 



3ft 2 W 
&, (70) 



(71) 



To see that there are only two characteristic variables 
here one must contract with the vector Sj. Finally the 
remaining characteristic variables are 

± ^ft^^ - ±^r, lm b^]%m (72) 

where [ ] denotes that the object is trace- free in down- 
stairs indices. For the highest frequency mode the Dq 
operator in the principal symbol vanishes. However, the 
symbol still has a complete set of characteristic variables. 
The conserved quantity may then be constructed as be- 
fore with a sum over the grid-modes. The conserved 
quantity is obviously a norm since it contains every grid- 
mode, and is therefore equivalent to the standard 

IMlLo+.pg = WiijWl + \\ a \\l + IIAIIfc + 



wml 



£(IP+*7tfl 



\D +k a 



k =i 



(73) 



in physical space. 

As the calculation does not rely in any significant way 
on the flat background, it should be simple to extend 
to the case in which one linearizes around an arbitrary, 
constant in space background. It may then be possible 
to extend to the case with variable coefficients in space 
following (6(| . We also anticipate no problems in extend- 
ing the calculation to higher order finite difference (FD) 
approximations. In our numerical tests in Sec. II VI we do 
not perform evolutions with this discretization. 



D. The algebraic constraints 

In this section we demonstrate that the numerical sta- 
bility of standard numerical implementation of the lin- 
earized LaSh (and BSSN) systems, which includes the 
conformal decomposition of the evolved variables, de- 
pends only upon the analysis of the previous section. In 
order to do so, we show that there is a one-one correspon- 
dence between solutions of the original and decomposed 
systems. 

In the linear regime the conformal decomposition is 
simply a linear combination of the undecomposed vari- 
ables subject to linear constraints. Consider the semi- 
discrete system under such a decomposition. Start by 
defining the decomposed state vector on a time slice by 
v = Tu. Here and in what follows we suppress spatial 
indices. Assume that u has m elements. Then T is an 
I x m matrix, v has I elements. We denote the pseudo- 
inverse of T by S, a matrix which maps from the image 
of T in R l back to R m such that 



ST = I n 



(74) 



If the evolution equations for u are given by Pu and 
those for the decomposed variables are Pv then the two 
are related as P = TPS. Denote by _L the projection 
operator which maps to the m dimensional hypersurface 
in R/ on which the algebraic constraints are satisfied. 
The algebraic constraints are 



C = v- _L v. 



(75) 



Consider first the semi-discrete system. Suppose that 
at a given time the constraints are satisfied. Then we 
find 



d t C = Pv- ±Pv = TPu- _L TPu = 0. 



(76) 



where the last equality holds because directly after the 
application of T the algebraic constraints are satis- 
fied, and therefore the projection operator does nothing. 
Therefore in the semi-discrete system if the constraints 
are satisfied initially they remain so, and there is a one- 
one correspondence between solutions of the decomposed 
and undecomposed systems. 

For the fully discrete system we take an explicit poly- 
nomial time- integrator Q for the undecomposed variables 
and the modified time integrator _L Q for the decomposed 
system. Now consider the difference between constraint 
satisfying data u n and v n = Tu n integrated with the two 
methods. For brevity we subsume the timestep At into 
P and P. One finds that 

v n+l _ Tu n+1 =± Q[p][ w »»] _ TQ[P][u n ] 

=_L Q[TPS][Tu n ] - TQ[P][u n ] 

= L TQ[P][u n ]-TQ[P][u n ] 

= TQ{P}[u n }-TQ[P]{u n }=0, (77) 

where we have used linearity of the system, polynomiality 
of the time-integrator and the fact that directly after the 
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application of T the algebraic constraints are automati- 
cally satisfied, so the projection operator does nothing. 
Thus the two integration methods are equivalent as de- 
sired. Note that in these calculations the modified time 
integrator _L Q could be replaced by Q since the unpro- 
jected time-step introduces no constraint violation. We 
have verified these calculations by explicitly comparing 
evolutions of the linearized conformal LaSh system with 
and without constraint projection. We prefer to discuss 
the natural linearization of the non-linear method, which 
includes the projection. In the non-linear case the unpro- 
jected timestep can introduce algebraic constraint viola- 
tions. 



IV. NUMERICAL EXPERIMENTS 

The LaSh system is implemented inside the Lean 
code [57f which is based on the Cactus computational 
toolkit [67j and the mesh refinement package Carpet 
[68l [o9| . Initial data is constructed by solving the 
constraint equations with the TwoPunctures spectral 
solver provided by [70l ]. 

We perform the following set of numerical evolutions: 

Robust stability tests We perform a subset of the so- 
called apples with apples tests [7l[ to demonstrate 
numerically that the evolution system is formally 
stable with various choices of the densitization pa- 
rameters; 

Puncture stability We evolve a single BH with differ- 
ent choices of the densitization parameters to estab- 
lish what restriction is placed on them by insisting 
on long-term stable puncture simulations; 

Head-on collisions We compare BSSN evolutions of 
the head on collision of two BHs with those per- 
formed with LaSh. We focus on consistency of the 
extracted physics at finite resolution; 

Binary black hole inspiral We compare BSSN evolu- 
tions of inspiraling BHs (Goddard RI [72]) with 
those performed with LaSh. We consider the com- 
putational costs of simulations with the two sys- 
tems as well as consistency of the results. 



A. Robust stability 

In Sec. IIII CI we did not succeed in demonstrating 
numerical stability of the LaSh system with the punc- 
ture gauge using standard discretization. We there- 
fore perform robust stability tests following the method 
in Ref. [Til ] . The numerical domain is given by 
-0.5 < x < 0.5, -0.06/A < y < 0.06/A and 
— 0.06/A < z < 0.06/A with periodic boundary con- 
ditions. We use three different resolutions h = 0.02/A, 
where A c = 1, A m = 2 and A/ = 4. The expansion in 



the y- and z- direction incorporates the three grid points 
required for fourth order FD stencils. The initial data are 
given by small perturbations of the Minkowski spacetime 

Jij = Vij + e»j • (78) 

The tij are independent random numbers in the range 
(-10" 10 /A,10" 10 /A), so that terms of the order 0(e 2 ) 
are below round-off accuracy. This means that the evo- 
lution remains in the linear regime unless instabilities 
occur. 

We monitor the performance of each simulation by 
calculating the maximum norm of the Hamiltonian con- 
straint as a function of time. For this study we focus on 
three choices of the densitization parameters (jiq,tik) = 
{(0, 0), (0.5, -0.5), (-0.5, 0.5)}. The results of the robust 
stability test are plotted in Fig. [TJ For all choices of 
i n Q, n K), including the BSSN scaling (0,0), we obtain 
stable evolutions. 



B. Puncture stability 

We next perform evolutions of a single puncture, 
studying a wide range of non-trivial densitization pa- 
rameters. The hyperbolicity analysis of the continuum 
LaSh scheme presented in Sec. IIII CI is not affected by 
the choice of densitization parameters provided that the 
algebraic constraints are enforced. In the previous sec- 
tion we have seen that various choices of the densitization 
parameters yield evolution systems that are numerically 
stable. Here we demonstrate that those parameters must 
be chosen more carefully in order to achieve long-term 
evolutions of puncture data. We evolve a single, non- 
rotating BH until t = 500 M. The BH is initially given 
by two punctures with mass parameter mi,2 = 0.5M lo- 
cated at z = ±10 _5 M. Using the notation of Sec. II E 
of Ref. [13] , the grid setup is given in units of M by 

{(96, 48, 24, 12, 6, 2, 1, 0.5), 1/32} . 

We vary both tiq and uk in the interval [—1,1] in steps 
of An = 0.1. The lifetimes T; of the simulations are 
determined as functions of the densitization parameters. 
The first occurrence of "nans" in the right-hand side of 
the densitized lapse Q is used as a measure of the lifetime 
whenever the simulations did not survive for the entire 
evolution. In Fig. [5] we show the results of this parameter 
study as a contour plot. In particular, a single puncture 
can be evolved for at least t = 500 M using the LaSh 
system with parameters indicated by the light blue area 
in the figure. Negative values of the lapse densitization 
parameter uq < —0.3 combined with positive values of 
the curvature densitization parameter tik > let the 
simulations crash after a short time. In contrast, long 
term stable evolutions are obtained for the parameter 
range n Q <G [-0.3, 0.9], n K e [-1, 0], including the BSSN 
scaling uq — uk — 0. 

We can partially understand this behaviour by consid- 
ering single puncture initial data and their influence on 
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FIG. 1: (Color online) Apples with apples stability test using a low (50 grid points, black solid lines), medium (100 grid 
points, red dashed lines) and high (200 grid points, green dashed-dotted lines) resolution. The tests were performed with the 
densitization parameters (nQ,n,K) = {(0,0), (0.5, —0.5), (—0.5,0.5)} from left to right, respectively. 
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FIG. 2: (Color online) Contourplot of the lifetime T; as func- 
tion of the densitization parameter uq and tik- Areas colored 
in dark blue indicate a short lifetime whereas light blue col- 
oring stands for a lifetime of at least Ti = 500A/. 



tial data they reduce to 

9t7«=0, dtX = 0, (80) 

d t A l 3 =x^ nK/2 [D t D ] u + atfjf 1 * , (81) 

d t K = - X ~ 3nK/2 D l D ia , (82) 

d t T 1 =0 . (83) 

Insisting on initially regular evolved variables at the 
puncture, Eqs. (|80M83p require tik < 0, also in agree- 
ment with our study; numerical experiments violating 
this condition immediately fail. 

For further illustration we plot in Fig. [3] the time 
derivatives of the densitized lapse Q and the trace of 
the extrinsic curvature K after an evolution time of 
t = 100 M. As we simoultaneously increase uq and 
decrease tik, we obtain smoother profiles. Note that the 
BSSN case uq = — tik produces the steepest gradients 
in this comparison. A systematic study of the excep- 
tionally benign behaviour of a non-trivial densitization 
on the accuracy of 10-15 orbit simulations, especially of 
spinning, precessing binaries, is beyond the scope of this 
paper. Our results may, however, point at fertile ground 
for future research of the LaSh system. 



the evolution equations (|16M20[) . On the initial slice the 
densitized lapse Q is given by 



Qo = X^o = X |( " Q+1/3) , (79) 



corresponding to a pre collapsed lapse a. Since x vanishes 
at the puncture we require nq > - j to obtain a regular 
densitized lapse Qq on the initial timeslice, in agreement 
with the findings of our parameter study; simulations 
with uq < —0.3 crash immediately. Next consider the 
evolution equations on the initial timeslice. For our ini- 



C. Head-On Collisions 

In this section we study in depth the stability proper- 
ties of numerical simulations of equal-mass head-on col- 
lisions performed with the LaSh system. For this pur- 
pose we evolve model BL2 of Table II in [57| . i.e. two 
non-spinning holes with irreducible mass M- lvV: i = 0.5 M 
starting from rest at Zip = ±5.12 M. The computational 
grid consists of a set of nested refinement levels given in 
units of M by 

{(256, 128, 72, 32, 16) x (4, 2, 1), h} , 

where we have usually chosen h — M/48, unless denoted 
otherwise. We consider the densitization parameters 
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FIG. 3: (Color online) Right hand sides of the densitized lapse Q (left panel) and of the trace of the extrin- 
sic curvature K (right panel) after an evolution time of t = 100M. We take parameter pairs (nQ,riK) = 
{(0, 0), (0.2, -0.2), (0.4, -0.4), (0.6, -0.6), (0.8, -0.8), (1, -1)}. 



Run 


Grid Setup 


d/M 


ng 


n K 


W 4 E rad /M 


HDl c 


{(256,128,72,32,16) x (4,2,1) 


h = 


1/40} 


10.24 


0.0 


0.0 


5.51 


HDl m 


{(256,128,72,32,16) x (4,2,1) 


h = 


1/44} 


10.24 


0.0 


0.0 


5.52 


HDl f 


{(256,128,72,32,16) x (4,2,1) 


h = 


1/48} 


10.24 


0.0 


0.0 


5.53 


HD2 


{(256,128,72,32,16) x (4,2,1) 


h = 


1/48} 


10.24 


0.2 


-0.2 


5.53 


HD3 C 


{(256,128,72,32,16) x (4,2,1) 


h = 


1/40} 


10.24 


0.4 


-0.4 


5.51 


HD3 m 


{(256,128,72,32,16) x (4,2,1) 


h = 


1/44} 


10.24 


0.4 


-0.4 


5.52 


HD3f 


{(256,128,72,32,16) x (4,2,1) 


h = 


1/48} 


10.24 


0.4 


-0.4 


5.53 


HD4 


{(256,128,72,32,16) x (4,2,1) 


h = 


1/48} 


10.24 


0.6 


-0.6 


5.53 



TABLE I: Grid structure and physical initial parameters of the simulations of a head-on collsion of an equal mass BH binary. 
The grid setup is given in terms of the radii of the individual refinement levels as well as the resolution near the punctures h 
(see Sec. II E in [571 ] for details). The table further shows the initial coordinate separation d/M of the two punctures. E ra d/M 
is the fraction of the total BH mass that is radiated as gravitational waves. All parameters are given in units of the total BH 
mass M = Mi + M 2 . 



(n Q ,n K ) = {(0,0), (0.2, -0.2), (0.4, -0.4), (0.6, -0.6)}, 
denoted as models HD1 - HD4 in Table HI For mod- 
els HD1 - HD3 we have chosen the T-driver shift con- 
ditions ([50)1 with (/is,£i,£2)?7) = (1, 0, 0, 1), whereas in 
case of model H DA the T-driver shift conditions (|30| have 
been taken with (jjl s , £i, £2) v) — (3/4,1,1,1) [73. In- 
formation about gravitational waves emitted during the 
plunge has been obtained by the Newman-Penrose scalar 
^4. In Fig. [4] we present the real part of the dominant 
mode , J , 20j rescaled by the extraction radius r ex — 60 M, 
for models HD1-HDA. Note, that the imaginary part of 
^4 vanishes due to symmetry. We find that the wave- 
forms generated by the different models agree well. Wc 
study the convergence of models HD1 and HD3 by us- 
ing three different resolutions h c = M/40, h m = M/44 
and hf = A//48 referred to as coarse, medium and high 
resolution. The differences of the t — 2, m — mode 
of the resulting gravitational radiation are displayed in 
Fig. [S]and demonstrate overall fourth order convergence 
for both models. We estimate the discretization error at 
high resolution in the waveforms ^20 to be 0.4%, similar 



to the error reported in [57[ for the corresponding BSSN 
evolutions. 

The amount of energy that is radiated throughout 
the head-on collision computed from, e.g., Eq. (22) in 
Ref. HH (see also 0) is E rad /M = 0.0553% for models 
HDlf, HD2, HDBf and HDA, again in excellent agree- 
ment with Ref. [57| ■ We estimate the discretization error 
in the radiated energy to be 0.4% and the error due to 
finite extraction radius to be 1.6%. 

As for single BH evolutions, we observe smoother time 
derivatives of Q and K for non-vanishing choices of uq 
and riK ■ We illustrate this behaviour in Fig. |5] which 
shows the time derivatives along the z axis obtained for 
different values of uq and uk at t = 10 M. 

D. Inspiraling Black-Holes 

In this section we will demonstrate how BBHs can be 
evolved successfully using the LaSh formulation of the 
3+1 Einstein equation in combination with the moving 
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Run 


hj volut ion 
Scheme 


Grid Setup 




UK 


10 2 E rad /M 


BSSN 


BSSN 


{(256, 128, 64, 24, 12, 6) x (1.5, 0.75), 1/48} 


- 


- 


- 


LaSho 


LaSh 


{(256, 128, 64, 24, 12, 6) x (1.5, 0.75), 1/48} 


0.4 


-0.4 




LaSh c 


LaSh 


{(256, 128, 64, 24, 12, 6) x (1.5, 0.75), 1/52} 


0.4 


-0.4 


3.69 


LaSh m 


LaSh 


{(256, 128, 64, 24, 12, 6) x (1.5, 0.75), 1/56} 


0.4 


-0.4 


3.68 


LaShf 


LaSh 


{(256, 128, 64, 24, 12, 6) x (1.5, 0.75), 1/60} 


0.4 


-0.4 


3.67 



TABLE II: Grid structure, evolution system and initial parameters of the simulations of quasi-circular inspirals. The grid setup 
is given in terms of the radii in units of M of the individual refinement levels as well as the resolution near the punctures h 
(see Sec. II E in [57]] for details). In case of the LaSh scheme we also specify the densitization parameters nq and uk- The 
final column lists the radiated energy 15 ra d extracted at r cx = 60 M for models LaSh c -LaShf . Models BSSN and LaSho have 
only been run until t = 50M in order to compare their computational cost. 
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FIG. 4: (Color online) Real part of r ox M\['20, the dimension- 
less Newman-Penrose scalar, where r ex — 60M, for model 
HDlf (black solid line), HD2 (red dashed line), HD3f (green 
dashed-dotted line) and HD4 (blue dotted line). 



puncture approach. For this purpose we consider the ini- 
tial configuration labeled Rl in Table I of Ref. 72] . This 
configuration represents a non-spinning, equal-mass bi- 
nary with a total ADM mass of M — 0.9957 in code units. 
The bare-mass parameters are mi^ = 0.483 and the BHs 
start at position x±.2 = ±3.257 with linear momentum 
Pi 2 = ±0.133 in the y-direction. The specifications of 
the grid setup, in the notation of Sec. II E of Ref. [57j . 
are given m Table UU For this model we have used 
the T-driver shift condition (|30l) w ith (/j, s , £i, £2, v) — 
(3/4, 1, 1, 1) as suggested in Ref. [73|. As before, we study 
the convergence properties by performing simulations of 
model LaSh with resolutions h c = 1/52, h m — 1/56 and 
hf = 1/60. In the left panel of Fig. [7] we present the 
real part of the I — 2, m = 2 mode of ^4, extracted 
at r ex — 40M, obtained by models LaSh c , LaSh m and 
LaSh c . The right panel of Fig.[7]shows the differences be- 
tween the coarse and medium and medium and high res- 
olutions of the amplitude (upper panel) and phase (bot- 
tom panel). The latter differences have been rescaled by 
the factor Q4 = 1.43 corresponding to fourth order con- 



vergence. The resulting discretization error in amplitude 
and phase are AA/A < 1% and Acf> < 0.1 rad. 

The energy radiated in gravitational waves is 
E rad /M = 3.67 ± 0.13% for the high resolution run 
LaSh f in Table |TT] which is in good agreement with the 
BSSN results of Ref. [E3. 



Run 


mem. 

[GByte] 


tr 

[CPUhours] 


V 

[M/hour] 


BSSN 


55 


290 


4.2 


LaSho 


70 


430 


2.9 


modLaSh 


55 


335 


3.7 



TABLE III: The required memory mem., the total runtime t r 
in CPUhours and the average speed v in units of physical time 
M per real time hour of the test simulations using the BSSN 
(model BSSN in Table OH), the original LaSh (model LaSho 
in Table [TTJ> and the modified LaSh scheme. The simulations 
have been run for t = 50A/ using 24 processors. 



Finally, we compare the computational performance of 
both, the BSSN and LaSh evolution scheme. For this pur- 
pose we have evolved models LaSho and BSSN until t = 
50M on the Magerit cluster [74( in Madrid which is part 
of the Spanish Supercomputing Network [75]. Magerit 
uses PowerPC-970FC processors running at 2.2GHz. 
The required memory, runtime and average speed ob- 
tained for 24 processors are shown in Table [TTTl The orig- 
inal LaSh system requires about 30 % more memory than 
the BSSN system and is about a factor 1.4 slower. The 
overhead of the LaSh system is not unexpected. First, 
the LaSh system involves a larger number of grid func- 
tions; the tracefree part of the extrinsic curvature A 1 j 
is not symmetric and thus requires 9 independent com- 
ponents instead of 6 for the BSSN variable Ay . Second, 
the densitization of variables requires extra variables and 
involved more complicated expressions on the right hand 
sides of the corresponding evolution equations. These ef- 
fects can be partly eliminated, however, without loosing 
the appealing properties of the LaSh system. For this 
purpose, we have tested a modified version of the origi- 
nal LaSh system, denoted as modLaSh in Table HTT1 Here 
we evolve Ay instead of the trace-free part of the extrin- 
sic curvature with mixed indices A*. As expected, this 
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FIG. 5: (Color online) Convergence analysis of the real part of ^20 of the Newman Penrose scalar $4, re-scaled by the extraction 
radius r ex = 40M, for models HD1 (left panel) and HD3 (right panel) in Table U We show the difference between the low 
and medium resolution (black solid line) and the medium and high resolution (red dashed line). The latter has been amplified 
by a factor of Q — 1.58 expected for fourth order convergence 




FIG. 6: Right hand sides of the densitized lapse Q (left panel) and of the trace of the extrinsic curvature K (right panel) after 
an evolution time of t = 10M for models HDlf (solid line), HD2 (dashed line), HD3f (dashed-dotted line) and HD4 (dotted 
line). 



modification equals the BSSN system in memory require- 
ments and also significantly reduces the computational 
costs relative to the original LaSh system. At the same 
time, however, modLaSh preserves the flexibility that 
has enabled us to obtain smoother behaviour of the vari- 
ables close to the puncture as compared with the BSSN 
scheme. 



V. CONCLUSIONS 

Motivated by a desire to better understand which 
are the important ingredients of the moving puncture 
method, we have studied the LaSh formulation of the 



Einstein equations. Provided that the algebraic con- 
straints of the system are imposed the formulation is 
equivalent to BSSN. Therefore we have investigated how 
the choice of evolved variables effects the success of nu- 
merical simulations of puncture initial data. The change 
of variable is parametrized by the densitization parame- 
ters (uq, n K ). 

We started by demonstrating that LaSh is formally 
numerically stable when linearized around flat space for 
arbitrary densitization parameters, with fixed shift and 
densitized lapse. A special case of this calculation is the 
numerical stability of BSSN. We attempted to show nu- 
merical stability of the system coupled to the puncture 
gauge, but find that the required calculations are too 
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FIG. 7: (Color online) Left: Real part of the dominant I — 2, m = 2 mode of the dimensionless Newman-Penrose scalar 
rM3i(\I r 4) ) where the extraction radius is r ex — 40M. The waveforms are shown for models LaSh c (black solid line), LaSh m 
(red dashed line) and LaShf (green dash-dotted line). Right: Convergence analysis of the Amplitude (upper panel) and phase 
(bottom panel) of the dominant I = 2, m = 2 mode of the Newman Penrose scalar We show the differences between the 
coarse and medium resolution (black solid line) and medium and high resolution (red dashed line). The latter difference has 
been amplified by Q4 — 1.43, indicating fourth order convergence. 



complicated even for computer algebra unless we move 
away from the standard discretization. 

We performed four types of numerical tests. The first 
class of tests includes robust stability test, specifically the 
so-called apples with apples tests. We find that the LaSh 
formulation is numerically stable for various choices of 
the densitization parameters. Next, we found that long 
term stable evolutions of single BH spacetimes require a 
more careful choice of the densitization parameters. It is 
interesting to note that the parameter choice correspond- 
ing to the BSSN system is located near the edge of the 
permissible range. Moreover, we have identified parame- 
ter choices which result in smoother profiles of the time 
derivatives of the evolution variables near the puncture 
as compared with the BSSN case. It will be interesting to 
investigate the impact of this behaviour on the accuracy 
of inspiral simulations lasting 10 — 15 orbits. While such 
a study is beyond the scope of this paper, it may pro- 
vide fertile ground for direct application of the results 
presented in this work. Furthermore, preliminary tests 
of higher dimensional BH spacetime evolutions indicate 
that the generalized BSSN formulation helps overcoming 
stability problems that have been encountered in D > 6 
[49l | . A more detailed analysis of this application will be 
presented elsewhere [7t| . 

We have further evolved head-on collisions as well as 
quasi-circular inspirals of binary BHs. For both cases, 
we have achieved long-term stable evolutions for a wide 
range of non-trivial parameter choices (tiq, Uk)- The 
evolutions produce convergent waveforms consistent with 
the BSSN results and comparable accuracy. As men- 
tioned above, we plan to compare the accuracy of both 
systems for more demanding inspiral simulations in fu- 



ture work. In any case, the binary simulations confirm 
the above finding that non vanishing values of uq and uk 
facilitate evolutions with smoother profiles of the evolu- 
tion variables in the neighborhood of each puncture. 

In summary, our results highlight the importance of 
the choice of variables for numerical calculations aside 
from any continuum PDE considerations. This opens up 
the possibility of significantly reducing errors in simula- 
tions of astrophysical binaries with large spins or mass ra- 
tios and also overcome stability issues reported for higher 
dimensional BH simulations 491. 
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